Design of ion channel blocking, toxin-like Kunitz inhibitor peptides from the tapeworm, Echinococcus granulosus, with potential anti-cancer activity

Over-expression of K+ channels has been reported in human cancers and is associated with the poor prognosis of several malignancies. EAG1, a particular potassium ion channel, is widely expressed in the brain but poorly expressed in other normal tissues. Kunitz proteins are dominant in metazoan including the dog tapeworm, Echinococcus granulosus. Using computational analyses on one A-type potassium channel, EAG1, and in vitro cellular methods, including major cancer cell biomarkers expression, immunocytochemistry and whole-cell patch clamp, we demonstrated the anti-tumor activity of three synthetic small peptides derived from E. granulosus Kunitz4 protease inhibitors. Experiments showed induced significant apoptosis and inhibition of proliferation in both cancer cell lines via disruption in cell-cycle transition from the G0/G1 to S phase. Western blotting showed that the levels of cell cycle-related proteins including P27 and P53 were altered upon kunitz4-a and kunitz4-c treatment. Patch clamp analysis demonstrated a significant increase in spontaneous firing frequency in Purkinje neurons, and exposure to kunitz4-c was associated with an increase in the number of rebound action potentials after hyperpolarized current. This noteworthy component in nature could act as an ion channel blocker and is a potential candidate for cancer chemotherapy based on potassium channel blockage.


Results
Sequence analysis of EgKU4 protein and model validation. The current study predicted the 3D structure of Kunitz4 using the full-length mature sequence of E. granulosus Kunitz4, the constituent of 84 amino acids, and the sequence of 1dtx, 1dem, 1dtk, and1tfx ( Figure S1). Modeller 9.22 was employed to compute structural models using the BLOSUM62 similarity matrix. The top 10 from among 10,000 constructed models, based on the lowest discrete optimized potential energy (DOPE) score, were selected for validation of stereochemical features. Valuation of model reliability was done in terms of Z-scores. Analysis of Ramachandran plot of the 3D structures from RAMPAGE showed that 94.7%, 3%, 2%, and 0% of residues were placed in the most favored, additionally allowed, generously allowed, and disallowed regions, respectively (Fig. S2, Supplementary Table 1). The z-score for the predicted model was − 4.4 with the ProSA web server, indicating an acceptable model quality compared with experimentally validated protein structures ( Figure S2e). Additionally, the results of a secondary structure alignment demonstrated a significant secondary structure compatibility between Kunitz4 and the template model ( Figure S2c). To assess the environment of each residue in the model, the compatibility of 3D-1D structures was evaluated using Verify 3D score.
Comparable interaction of Kunitz4 and EAG1 in three distinct sites. Findings of the study indicated that the putative binding site of EAG1, located close to the selectivity filter between S5 and S6 segments and consisting of Tyr423, Met478, Tyr347, Lys340, Glu346, Ala350, Thr430, His364, Tyr344, and His343, belongs to the extracellular opening of the pore (Fig. 2d). The docking poses revealed three distinct binding sites interacting with Kunitz4 on the N-terminal and C-terminal of the EAG1 protein comprised of ranges of 18-35, 46-55, and 65-84. These sites were probably involved in the channel interaction, so these amino acids were selected for the design of peptides (Fig. 2a-c). Regarding the binding mode, analysis of the docking revealed that the Kunitz4 binding was centered, the number of hydrogen bonds was considerable, and the hydrophobic interactions were more dominant than other interactions. Almost all of the interaction patterns showed that they shared several common conserved residues based on multiple sequence alignment in their interactions. Based on molecular docking, consistency was found between the predicted active and passive interfacing amino acids and the functional amino acids in either EAG1 or Kunitz4. The lowest energy of binding for the EAG1-Kunitz4 complex was − 3149. 6 Fig. 2e places Kunitz4 Arg37 and Gly62 side chains in close proximity to Asn441 (A chain) of EAG1; with a fully elongated Kunitz4, Asn23 formed a hydrogen bond to Tyr344 (D chain) of EAG1 with a length of 3.34 Å. Dramatically, the peptides revealed a very low binding free energy in the complex with EAG1 (Table 1). Peptide-EAG1 docking turned out to be compatible with the mode and affinity of the Kunitz4 and EAG1 interaction ( Figure S3).

Stability of predicted model and Kunitz4-EAG1 complex after MD simulation.
The structural stability of the model was studied using molecular dynamics (MD) simulation. The output was analyzed to determine conformational changes against the initial structure by root mean square deviation (RMSD) based on the structure of the backbone. It is used in protein structure predictions to assess the match between the modeled and X-ray/NMR protein structure and often computed using only the Cα or backbone heavy atoms of each amino acid. The RMSD plot of the Kunitz4 structure fluctuated around 28,000 ps of the simulation and then reached a steady state within 0.3-0.4 nm of RMSD up to 50,000 ps (Fig. 3a). Radius of gyration (Rg) was calculated to estimate the compression and folding of the protein structure during the simulation (Fig. 3b). The Rg value of Kunitz4 increased slightly from 1.30 nm at the beginning to 1.36 nm around the 26,000th ps of simulation, and stability was maintained until the end of the simulation. This demonstrated that the Kunitz4 structure reached a stable and compact conformation during the simulation. To recognize the dynamics and function relationships of protein movements, the residue-based root mean square fluctuation (RMSF) values of the backbone Cα were compared. A large RMSF value indicated a flexible region with motion, while a low RMSF value proposed a rigid region and minimal motion. During simulation, the RMSF of Kunitz4 was less than 0.2 nm for all residues, indicating low changes in the protein structure. Molecular dynamics simulations for Kunitz4-EAG1 indicated the stability of the formed intermolecular interactions during the MD simulation, RMSD showing the stability of the complex during simulation with fluctuation of 0.3-0.8 nm (Fig. 3c). The Rg plot showed that Kunitz4-EAG1 became compact and relaxed at 30,000th ps of the simulation (Fig. 3d). The RMSF plot of EAG1 revealed slight fluctuations during the simulation. This could be interpreted as the higher stability of the four chains compared to Kunitz4, indicating approximately the same fluctuations in both EAG1 binding sites (Fig. 3f). In contrast, the RMSF plot of Kunitz4 in complex with EAG1 www.nature.com/scientificreports/ showed that the selected region for peptides fluctuated slightly in complex in comparison with free Kunitz4, at least in the two binding sites of EAG1 and Kunitz4, while binding sites between 64 and 75 exhibited minimal fluctuations (Fig. 3e). The number of hydrogen bonds is an index of the quality of the interaction. The total number of hydrogen bonds in the binding site during simulation varied between 40 and 60, and these findings validated the previous docking results (Fig. 4a). Accessible surface area solvent (SASA) was calculated for the interface of EAG1 and Kunitz4. Higher SASA values indicated more exposed amino acids in solvent, whereas lower values indicated more concealed residues in interaction. SASA plot analysis determined that the accessible area at the start of simulation was around 285 nm 2 , while the area was ultimately lowered to the range of 260-265 nm 2 at the end of simulation (Fig. 4b).
Molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA) calculation and binding energies. The EAG1-Kunitz4 complex was analyzed by considering residues contributed favorably to ΔE binding by means of energy decomposition calculation. Hot spots were identified by looking at residues that contributed significantly to complex stabilization. The area that showed the highest stabilization was identified in three specific regions located in the contact surface between EAG1 and Kunitz4. Supplementary Dataset file 2 shows residues of Kunitz4 that contributed to the stabilization of the complex. As shown in Fig. 5, in fragment kunitz4-a, the residues Arg21, Val22, Asn26, Leu27, Pro28, Ile29, Lys31, Gln33, and Arg35, in fragment kunitz4-b the residues Pro45, Ser46, Lys47, Arg52, and Phe53, and in fragment kunitz4-c, the residues Arg64, Lys66, Lys68, Arg69, Lys72, and Arg73, had more than 10 kj/mol of free energy contributed to their binding to the protein.  were evaluated based on their binding mode to EAG1 by ClusPro and HADDOCK web servers. Models with interactions were selected based on the lowest binding energy, higher hydrogen bond, and higher binding affinity to EAG1. Complexes formed from kunitz4-a and kunitz4-c were more stable with lower binding energy values of ΔG = − 8.1 kcal/mol and ΔG = − 7.4 kcal/mol, respectively, compared to the kunitz4-b complex with ΔG = − 6.8 kcal/mol. Likewise, measuring dissociation constant (Kd) showed that peptides a and c had lower kd compared to peptide a, with 1.5 × 10 -6 M and 2 × 10 -6 M vis-a-vis 6 × 10 -6 M, respectively (Table 1). These derived peptides also showed higher numbers of hydrogen bonds and lower binding energy compared with peptide b.
Peptides inhibit HepG2 and HT29 cell proliferation. MTT assays were used to test the effects of the three peptides on the proliferation of cancerous HT29 and HepG2 cell lines as well as the normal Hek293 cell  (Fig. 6). After 48 h, no significant suppression of normal Hek293 cells displayed in treatment with the peptide kunitz4-c up to 300 μM, however at 400 μM, between 25 and 35% decline in viability were obtained using all the peptides. The optimal amount of IC50 in cancer cells was obtained in the range equal to 251 µM for Hep G2 and 15.69 µM for HT29 in 48 h. Therefore, to investigate cytotoxicity effects of a 48-h study, comparing data of cancerous and normal cell lines was only considered at concentrations in the range of 0-400 μM (Fig. 6).
Immunocytochemistry staining. The expression of Ki-67 was analyzed in HT29, and HepG2 cell lines for kunitz4-a and kunitz4-c in 100 and 250 µM, respectively. The findings indicated a significant difference in the cells treated with small peptides compared to the untreated controls. As shown by densitometry analysis, the significantly reduced expression of Ki-67 was observed after treatment with the small peptides kunitz4-a and kunitz4-c compared to the untreated cells (p < 0.0001, Fig. 6-Panel B).
Kunitz4 peptides promote apoptotic cell death. To examine the apoptotic effects, HT29 and HepG2 cell lines were treated with Kunitz4 peptides and stained with annexin V and PI. As shown in Fig. 7, after 24 h of exposure to 100 μM of each peptide for HT29, and 250 μM for HepG2 and the normal Hek293 cells, a significant  www.nature.com/scientificreports/ increase in apoptosis was observed in HT29 (p = 0.005) and HepG2 (p = 0.001) cell lines. The ability of kunitz4c to promote apoptosis was significantly stronger on HepG2 (p = 0.0009) than on HT29 cells (p = 0.005). The peptides induced low apoptotic effects on the normal cell line compared to the cancerous cells, i.e. we reported 5.3% and 9.73% apoptosis for kunitz-b and kunitz-c, respectively, which are significantly lower than the corresponding values for HepG2 (19.25% and 50.92%) and Ht29 (22.08% and 16.46%). The results for kunitz-b and kunitz-c were more satisfactory than kunitz-a. The maximum amount of apoptosis level was observed in HepG2 occurred under kunitz4-c.
Derived peptides caused remarkable cell cycle arrest at the G0/G1 phase. To investigate cell cycle arrest involvement in the cell growth blockade, cell cycle alterations in HT29, HepG2, and Hek293 cells following Kunitz4 peptide treatments were evaluated. As shown in Fig. 7, all peptides caused a significant increase in both HT29 and HepG2 cell populations at the G0/G1 phase and S phase along with a decrease in cell population in G2/M compared to the control group (p = < 0.0001), while no remarkable cell cycle effects were found on normal Hek293 cells treated with peptides (p = 0.7217). These findings indicated that the small peptides inhibited the growth of HT29 and HepG2 cells with cell cycle arrest at the G0/G1 phase. Expression levels of the P27 (CDKN1B) gene in HepG2 and HT29 cell lines were increased 0.04 and 172 times for kunitz4-a and 6.5 and 188 times for kunitz4-c, respectively (Fig. 8a, b). However, following treatment with kunitz4-a, the expression level of CDK2 was found to be significantly decreased (p = 0.015) in the HT29 cell line Likewise, P53 expression, with an approximate molecular weight of 44 kDa, was markedly increased by kunitz4-a and kunitz4-c in HT29 and HepG2 cell lines, suggesting that both kunitz4-a and kunitz4-c could increase the intracellular expression level of the tumor suppressor P53 gene in tumor cells.
Action potentials were induced in Purkinje neurons using negative current pulses 520 ms in duration and − 0.5 to − 0.1 nA (0.1 nA increments) from V holding = − 60 mV. The first spike latency and rebound action potentials generated by negative current injections were measured as described previously. Purkinje cells from peptide-treated rats in the presence of a hyperpolarizing conditioning pulse (− 0.1 to − 0.3 nA) displayed a significant decrease in the first spike latency (Fig. 10). Exposure to the peptides was also associated with a significant increase in the number of rebound action potentials after the application of hyperpolarizing current pulse − 0.1 nA and a significant increase in the − 0.2 nA only for kunitz-c (Fig. 10). The voltage-gated K+ current in Purkinje cells was evoked with a step-up depolarization protocol (Fig. 11a, b). Briefly, the membrane potential was depolarized to + 30 mV (10 mV increments per 10 steps, duration = 100 ms) and subsequently restored to the original potential. Perfusion of 200 μM of each of the peptides (n = 6) produced a slightly significant reduction in peak amplitude ( Fig. 11c; command voltage = + 30: control = 3045 ± 68; kunitz4-b = 2685 ± 51.3, p = 0.012; kunitz4c = 2816 ± 43.1, p = 0.032) and the steady-state current (control = 2468 ± 39; kunitz4-a = 2115 ± 34.6, p = 0.045; kunitz4-c = 2031 ± 42.1, p = 0.04), whereas 200 μM of kunitz4-a (n = 6) produced a non-significant peak current.

Discussion
In the present study, we demonstrated the anti-tumor activity of Kunitz-type protease inhibitors derived from the tapeworm Echinococcus granulosus. Kunitz-type toxins occur naturally in many poisonous animals, such as snakes, sea anemones, and cone snails, and in the saliva of bloodsucking arthropods 22 . The potassium ion channel blocking role of Kunitz-type toxins has been well-documented in collinein-1, the venom of Crotalus durissus collilineatus, and Hirudin, derived from the leech, Hirudo medicinalis saliva 26,27 . The significance of potassium channels in tumor growth has been clarified in several studies on cancer biology 28 . Current evidence indicates that the upregulation of voltage-dependent potassium channels is connected to cancer hallmarks, and the ectopic expression of Kv10.1 in tumors occurs in the majority of human tumor types 12 . It has already been shown that the E. granulosus Kunitz family demonstrates peptidase inhibitory properties as well as channel blocking effects. Among different members of the family, EgKU-1 and EgKU-4 (Kunitz1 and Kunitz4) induce their effects through voltage-activated potassium channel blockade 24 . In the present study, we characterized the anti-cancer activity of synthetic peptides derived from E. granulosus Kunitz4 against two cancer cell lines using computational modeling as well as in vitro cellular methods and biological assays. In addition, protein-protein interactions and Kunitz4 affinity behavior for the EAG1 channel were investigated using molecular simulation. The current findings of modeling showed that synthetic peptides derived from Kunitz4 can induce significant effects equivalent to those of the whole protein.
As therapeutic agents for the incorporation of protein-protein interactions (PPIs), peptides are able to span a large contact surface area, making them suitable for binding to specific cell surface receptors with high specificity and potency; peptide toxins propose a better bioavailability and stability, such as G protein-coupled receptors or ion channels, and consequently, they trigger intracellular targets 29,30 . In designing the peptides, before preparing the complex, we first predicted and validated the proper 3D structure of Kunitz4 as well as the Kunitz4 binding site through multiple ways, such as 3D LigandSite, COACH, and CPORT servers.
Homology modeling of Kunitz4 was performed using 3D conformation of α-dendrotoxin in the venom of green mamba, Dendroaspis angusticeps (PDB code: 1dtx). High sequence homology between the toxin and Kunitz4 (e.g., 45% sequence identity between α -DTX and Kunitz4) indicates a similar main-chain fold. Kunitz4 was composed of two α helixes and two anti-parallel β sheets and random coils. Interestingly, sequence analysis of Kunitz4 and 1dtx indicated five cross-linked disulfide bonds between the conserved Cys residues, suggesting that cysteine residues are drastically important to maintaining the tertiary structure. The comparison of α-dendrotoxin and Kunitz4 showed that both molecules share amino acids that can be important for biological activity; for example, Leu27 and Lys47 in Kunitz4 are equivalent to Leu9 and Lys29 in α-dendrotoxin, respectively. These observations signify the potential of Kunitz4 as a cation-channel blocker 24 . Selected docking poses of Kunitz4 amino acids interacting with EAG1 indicated that these amino acids were mostly scattered over the Kunitz4 sequence. Three domains of Kunitz4 are important in interactions with the channel; amino acids within the 21-35 and 42-53 regions of Kunitz4 are involved in the interaction with EAG1, and the third area www.nature.com/scientificreports/ located between 64 and 76 from Kunitz4 is evolutionarily notable. Along with molecular dynamics simulation, trajectories analysis showed that the predicted structure of Kunitz4 was stable with low fluctuations during the simulation. Moreover, the RMSD plot of the complex of EAG1-Kunitz4 showed convergence after the initial fluctuations. Interestingly, the binding sites of Kunitz4 showed higher fluctuations in the free form compared to the complex with EAG1. Evaluation of the stability of the complex by molecular dynamic simulation showed that the designated peptides interacted effectively in the protein-channel interactions (Fig. 3). Different K channel-blocking toxins usually share common features essential for interaction with the channel, e.g., functional pairing of the combination of lysine and phenylalanine/tyrosine as hydrophobic residues. This allows the electrostatic interaction of the positive-charged lysine with the negative residues in the pore, causing the physical block of ion conduction, as was shown for Kunitz4 in the present study. Eleven amino acids from Kunitz4 made hydrogen bonds with Kv10.1, and 31 residues were involved in hydrophobic contacts with the channel; most important among them is the conserved lys47 binding to Met478 in the S5-S6 domain. The hydrogen bond triad formed by Arg21, Asn23, and Ile29 with Tyr 344 from the channel made a binding spot that could interrupt the regular octahedral geometry of the open pore.
The crystallography structure of voltage-dependent K+ channel blocking toxins indicated that the effective strategy required a highly conserved residue, Lys47, in the selectivity filter at the extracellular K + binding site 31 . Hotspot residues of EgKU4 computationally predicted with the MM-PBSA approach showed that Leu27 and Lys47 play fundamental roles in PPI (− 17 kJ/mol and -182 kJ/mol, respectively). Other critical residues constructing three peptides include Arg21, Lys31, Arg 35, Arg 52, Asn61, Arg 64, Lys66, Lys68, Arg 69, Lys72, and Arg 73. Considering the data, the three regions of Kunitz4 were selected for peptide design against the extracellular S5-S6 domain using hotspot residues prediction with MM-PBSA in GROMACS. As is the case in charybdotoxin isolated from Leiurus quinquestriatus, the scorpion toxin that fits into the pore of the channel in a lock and key posture 31 , peptides block EAG1 in the same manner ( Figure S3).
The first natural toxin inhibiting Kv10.1 (EAG1) was reported from scorpion toxin κ-hefutoxin 1, isolated from Heterometrus fulvipes 32 . The MIC50 value of the toxin against Candida albicans, Saccharomyces cerevisiae, and Fusarium oxysporum was found to be in the range of 18.8-37.7 µM. APETx4, a peptide derived from the sea anemone Anthopleura elegantissima, as a gating modifier, has been shown to induce inhibitory effects on EAG1 through binding to the S3-S4 region. Furthermore, it has displayed concentration-dependent proapoptotic and cytotoxic effects in LNCAP, SH-SY5Y, and MDA-MB-435S cancer cell lines 33 . Collinein-1, a snake venom thrombin-like enzymes (SVTLEs) reduced the viability of the human breast cancer cell line MCF7 but did not affect either the HepG2 cell line or MCF10A, the non-tumorigenic epithelial breast cell line. Arg79 residue is essential for the Kv inhibitory motif and interacts directly with the EAG1 channel selectivity filter 27 . In a xenograft model of liver cancer, Procyanidin B1, a natural composition from the grape seed, was identified as a specific inhibitor of Kv10.1; it inhibited tumor growth in vivo and suppressed proliferation and migration of hepatoma cells 34 . Herein, we evaluated the anti-proliferative, cytotoxic, and proapoptotic effects of Kunitz4-derived peptides on cancerous and normal cell lines. The peptides were shown to be able to induce a dose-dependent effect on both cancer cell lines; however, comparing to HepG2 and Ht29, very low anti-proliferative and cytotoxic effects were found on Hek293, the normal control cell line. It should be noted that the cytotoxic effect was remarkably higher on colorectal adenocarcinoma cell line HT29 than HepG2, the hepatocarcinoma cells, with lower IC50s. An increased apoptosis rate was also observed after treatment of the cells with certain concentrations of the peptides. As shown in Fig. 7, the cancer cells displayed typical apoptotic features; however, minimal changes were seen in noncancerous Hek293 cells.
Much evidence has shown that potassium channel activity is required for G1 progression of the cell cycle in different cell backgrounds, suggesting that potassium channel activity is necessary for early-stage cell proliferation through the G1 phase and facilitation of G2/M progression during the cell cycle. Moreover in cancer cells, the activity of EAG1 is high in the G1 phase and decreases when the cells enter the S phase 35 . The current findings showed a significant accumulation of cancer cells in the G0/G1 phase as well as decreases in cells in the G2/M and S phases in comparison to the control cells (Fig. 7). The results suggest that cell proliferation is suppressed by the peptides through cell cycle arrest at the G0/G1 phase, providing a chance to repair and discontinue proliferation of damaged cells.
Notably, high expression of P53 is involved in the regulation of cell cycle arrest. Potassium channels are involved in apoptosis by regulating cell volume; when the regulation fails, the outflow of K+ ions become too intense, water leaves the cells, and that results in a decrease in cell volume and subsequent cell death. P53 is a master regulator of DNA replication and cell death in cancer cells and controls the three checkpoints G1, S, and G2/M by regulating CDKs and cyclins as positive and P16, P27 and P21 as negative regulators of the cell cycle. In the present study, by comparing responses of the two cell lines with the corresponding controls, the synthetic peptides showed a significant impact on the cell cycle following upregulation of P53. The same trends were found in cell cycle arrest caused by the peptides. These results demonstrated that kunitz4-a and kunitz4-c could probably inhibit the downstream signaling through cell cycle arrest involving the up-regulation of P53 and P27. Our findings are in line with the fact that the tumor suppressor P53 gene induced cell cycle arrest at the G0/G1 phase. Therefore, p53-triggered inhibition of cell proliferation could be attributed to cell apoptosis and induction of cell cycle arrest 36 . The expression of P27 in both HepG2 and HT29 cells was remarkable in this study, and the distribution of HepG2 cells in the cell cycle showed that the peptides significantly increased the cells in the G0 and G1 phases, while cells in the S and G2-M phases were decreased. P27, a main cell cycle regulator, accelerates cell cycle development in HT29 and HepG2 cells; therefore, a 2-to 3-fold increase in P27 expression is adequate for obstructing the G1-S-phase and inhibiting cyclin-CDKs 37 . TEA, 4-AP, and glibenclamide blocking the Kv and K ATP channels in the U87 MG glioma cell line could inhibit tumor growth by arresting G0/G1 transition 38  www.nature.com/scientificreports/ evaluated by Ki-67 protein using ICC. The Ki-67 protein, with its relatively short half-life is present during all active phases of the cell cycle (G1, S, G2, and M). During anaphase and telophase, a firm decrease in Ki-67 levels occurs. Expression of the Ki-67 protein, as a marker of tumor aggressiveness, is associated with the proliferative activity of malignant tumor cell populations, as it has been shown in a number of studies on cancers of the cervix, lung, prostate, breast, soft tissues and the central nervous system 39,40 . The present study demonstrated that kunitz4-a and b decreased CDK4 mRNA in HepG2 cells, while only kunitz4-c caused an increase in the expression of P27; the other peptides had no significant effect on P27 in HepG2 cells (Fig. 8). In HT29, however, this condition was different, and all the peptides could produce significant amounts of P27 and reduce CDK2 expression with no considerable effect on CDK4. These findings indicate that the peptides inhibit cell proliferation by inducing G1 arrest through accumulative effects of P27 and P21, two well-known CDK inhibitors regulating cell proliferation41. Huang reported that treating the U251 glioma cell line with K ATP channel blockers resulted in cell cycle blockage in the G0/G1 phase 10 . In another study, G2/M cell cycle blockage and apoptotic cell death were documented following Ca2+ -activated Kv channel KCa3.1 blocking in GL261 glioma cells treated with temozolomide 42 . The present study provides evidence of change in gene expression at the RNA level of P53 and P27 caused by some peptides in tumor cells, however, this difference does not translate to the protein level which is expectable because of the extensive regulation of these proteins, especially P53, at the posttranslational level. It should also be noted that, other factors determining the specificity of the Kunitz4 peptides towards cancer cells need to be further studied. Cancer cells have more positive membrane potentials than normal cells. It has been demonstrated that cells are depolarized in the early G1 phase and hyperpolarized during the progress from the G1 into the S phase 41 . As shown in Figs. 7, 9, and 11, inhibition of K+ channels or cell blockade in the G1 phase is associated with membrane depolarization. Likewise, membrane depolarization initiated by an increase in the concentration of extracellular K+ simulates the effects of potassium channel blockers. Alterations in the firing characteristics of Purkinje cells are reflected in the functional properties of ion channels, which provide basic regulation of nerve excitability. The current results showed that treatment with Kunitz4 peptides increased the Purkinje neuronal intrinsic excitability, which is thought to be a function of inhibition potassium channels, especially by kunitz4-c. Purkinje neurons exposed to peptides exhibited significantly shorter duration of action potential and increased action potential frequency.
Treatment with peptides also induced a significant decrease in latency to the first spike that was associated with a significant increase in the rebound spike firing at the offset of hyperpolarization. This is in accordance with findings of previous studies which have revealed that the blockade of KV channels affects not only the first spike latency, but also the after hyperpolarization and frequency of spikes 43 . The shortening of the first spike latency could be attributed to the inhibition of the fast transient (A-type) channel, which is thought to play a key role in the timing of rebound action potentials 44 . The three peptides synthesized based on E. granulosus Kunitz4 showed anti-cancer activities against two human cancer cell lines, namely HepG2 and HT29. Cellular and molecular findings based on cell viability, cell cycle and apoptosis analyses, immunocytochemistry, major cancer cell gene involvement in cell progression, and whole-cell patch clamp recordings showed that the peptides have highly selective inhibitory effects on A type channels and are capable of suppressing tumor cell growth by inducing cell cycle arrest and apoptotic cell death. However specific actions of the peptides through the inhibition of K+ channels (specifically EAG1 channel) should be substantiated by the correlation of EAG1 protein level and activity with the effect of the peptides in several cell lines. The peptides characterized and evaluated in the present study have a promising potential for further in-depth cellular and molecular investigations in both in vitro and in vivo settings.

Methods
Sequence analysis of Kunitz4 protein and its homologs. The Kunitz4 protein sequence of Echinococcus granulosus G1 genotype was extracted from NCBI (APY21165.1), and the sequences of the dendrotoxin I from Dendroaspis polylepis polylepis (black mamba) (NCBI:P00979), α-dendrotoxin from the green mamba venom (NCBI:P00980), Dendroaspis angusticeps (eastern green mamba), (NCBI:P00980), homolog K from D. polylepis polylepis as dendrotoxin k (NCBI:P00981), and the Kunitz domain of tissue factor pathway inhibitor (NCBI:P10646) were used to compare Kunitz4 and sequence alignment analysis of the toxins as well as for the homology modeling of Echinococcus granulosus Kunitz4.
Kunitz4 structure determination. Based on the aforementioned templates, the three-dimensional protein structure of Kunitz4 was predicted by a model constructed from its primary amino acid sequence (NCBI: APY21165.1) using Modeller 9.22. Structure validation was performed using MolProbity 45 , ProSA 46 , ERRAT, and Vadar web servers 47 . The reliability of the homology model was assessed by calculating the Z-score and Ramachandran plot. The hydrophobicity plot was obtained in Discovery Studio 2.5 according to amino acid residue in the favored regions 48 .
Active site prediction. 3D LigandSite, CPORT (https:// milou. scien ce. uu. nl/ servi ces/ CPORT), and COACH (https:// zhang lab. ccmb. med. umich. edu/ COACH) were the programs utilized for predicting proteinprotein interface residues. Among the noted, six interface prediction methods were combined into a consensus predictor by CPORT, and the predictions were used as active and passive residues of Kunitz4 protein in molecular docking. www.nature.com/scientificreports/ (PDB ID: 5K7L) as a template to construct a continuous Kv10.1 channel structure. M-ZDOCK was utilized to predict the complementarity tetramer of Kv10.1 based on the structure of an unbound or partially bound monomer structure. Finally, each pose was weighted based on energy scores as described in detail by Chen and Weng. As a flexible docking approach, HADDOCK encodes information from predicted protein interfaces for the modeling of complexes 49 . We used Cluspro 2.0 in the rigid body docking step as another docking program which uses PIPER, an FFT-based docking algorithm. The Kunitz4 protein and EAG1 active site residues were defined for docking (Supplementary Dataset File 1). To this end, the highest scores were analyzed to identify the best channel-toxin-like conformation. Using Discovery Studio 2.5 48 and Ligplot+ 50 , all top-scoring conformations of every protein docking server were analyzed in terms of the number of hydrophobic and H-bonds. Finally, the top-ranking complex of EAG1-Kunitz4 models was subjected to 50,000 ps MD simulations and alanine mutation based on molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA). To check the molecular interaction of EAG1-peptides, molecular docking was carried out using ClusPro and HADDOCK. EAG1-peptide complexes were subsequently analyzed by Ligplot+ and Chimera. The complexes with the lowest interaction energy and the highest interaction bonds were selected for channel inhibition. Using the PRODIGY server, the binding affinity (ΔG) and dissociation constant (Kd) were calculated for selected peptide complexes.

Molecular dynamics simulations. The stability and conformation of the docked complexes and 3D
structure of Kunitz4 were determined using the program package GROMACS 5.1.4 51 along with gromos54a7 ff force field. The system was solvated with TIP3P water and simulated in an octahedron box with periodic boundary conditions (11 Å padding in each direction). Then the net charge system was neutralized by adding sodium and chloride ions to the protein. Minimization was performed with 5000 steps by the steepest integrator, and the maximum force was less than 1000 kJ/mol/nm. Subsequently, the system was equilibrated for 1 ns in the NVT and a constant pressure (NPT) of 1.0 bar. The LINCS algorithm was applied to all bonds involving hydrogen atoms during the heat and constant pressure steps. Force constants of 1000 kJ/mol/nm 2 were employed during the equilibration phase. Afterwards, the isothermal-isobaric MD simulation produced a 50,000-ps trajectory with a non-bonded cutoff distance of 14 Å and a leapfrog integrator with a 2.0 fs time step. All bonds involving hydrogen atoms were constrained by the LINCS algorithm during the heating and the constant pressure steps.
Using the Particle Mesh Ewald method, the electrostatic interaction energies were calculated, and the energies and coordinates of each atom were stored every 2 ps for subsequent analyses. Regarding the system behavior, RMSD, RMSF, Rg, SASA, and number of hydrogen bonds in the Kunitz4 and EAG1 interface were extracted from trajectories. The data was visualized using GraphPad Prism version 9.0.0.

Alanine mutation and binding free energy calculations.
To calculate free energy of binding (ΔG) and electrostatic energy, VdW energy, polar solvation energy, and surface accessible surface area, we used the Molecular Mechanics Poisson-Boltzmann Surface Area (G_MMPBSA) process executed in GROMACS 5.1.4 52 .
A total of 38,000 snapshots of structures at intervals of 10 ps were extracted from 50,000 frames of the MD runs and used to calculate the system enthalpy using MmPBSADecomp.py script provided along with the g_mmpbsa tool.
Peptide design. Three distinct sites from Kunitz4 that interact with EAG1 were selected for peptide design.
A peptide was constructed based on all possible amino acid substitutions. Key residues were recognized in Kunitz4-EAG1 interactions according to molecular docking studies, literature data mining, and multiple sequence alignment of the target to find conserved areas. Most processes involving peptide interactions were constructed using molecular dynamics-based methods as well as free energy calculations of snapshots equally spaced along a single dynamical trajectory 53 . Using the calculations, the experimental ΔΔG of binding with a mean error of ± 1 kcal/mol was predicted for the alanine mutations of hydrophobic and polar/charged residues with no buried salt bridges.
Evaluation of physicochemical properties of the peptides. The physiological/biochemical properties of the peptides, including molecular weight, isoelectric point, aliphatic index, water solubility, net charge at pH 7, grand average of hydropathicity (GRAVY), aggregation hot spots, and instability index, were calculated using ProtParam and PepCalc 54 . A negative GRAVY value indicates hydrophilicity, while a positive value indicates the hydrophobicity level of a protein. The best peptides with completely matched properties for in vitro experiments were selected based on the binding affinity energy, predicted dissociation constant (Kd), and stability (Table 2).  with 250 µM and HT29 cells with 100 µM of each peptide, total RNA was isolated using the TRIzol protocol. An equal amount of total RNA was used to synthesize cDNA using an Easy cDNA Synthesis Kit (Pars Tous, Mashhad) as instructed by the manufacturer. The expression of five major cancer cell biomarkers, namely P53, P27, P16, CDK2, and CDK4 genes, were determined by qRT-PCR. Two-step PCR conditions were used with the following thermal profile: 95 ℃ for 3 min, 40 cycles of 95 ℃ for 6 s, and 60 ℃ for 40 s. For each cell line, two internal standard genes were used as follows: GAPDH and PGK1 for HT29, GAPDH and TFG for HepG2, and HPRT and GAPDH for HEK293 55 (Supplementary Table 2). Data was statistically analyzed using the 2−∆∆Ct method 56 .
Immunocytochemistry staining. Immunocytochemical analysis was performed on the detached cells.
Smears were prepared on adhesive, positively-charged slides after 24 h of incubation of HepG2 and HT29 cells with 250 µM and 100 µM, respectively, of kunitz4-a and kunitz4-c. The slides were washed three times with PBS after 1 ml pure ice-cold acetone/methanol (1:1) was added to the slides at − 20 °C. Immunostaining was performed with 1:200 diluted mouse anti-Ki-67 monoclonal antibody (Zytomed, mouse anti-Ki-67, MSK018), and the slides were incubated for 30 min at room temperature. Subsequently, secondary HRP antibodies were used (Zytomed, Germany) with incubation times of 30 min and counterstained with hematoxylin (Dako, Waldenbronn, Germany) for a few seconds. The cells were analyzed and photographed at × 400 magnification using a digital microscope camera coupled with a computer (Olympus BX53; Olympus, Tokyo, Japan). To analyze immunoreactivity, three slides from each experimental group and 20 fields of each slide were examined to acquire the intensity of reactions.
Statistical analysis. The data analysis and plots were performed in GraphPad Prism version 9.0.0. The error bars on the figures correspond to the SEM. Statistical analysis was performed using the Student's t test or one-way ANOVA and two-way ANOVA method with Tukey's post-test, comparing all treatments to the negative controls and considering values of P values < 0.05 were considered significant. n represents the number of recordings.

Data availability
The used and/or analyzed datasets during the current study are available from the corresponding author on reasonable request.